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The granular silo is one of the many interesting illustrations of the thixotropic property of granular 
matter: a rapid flow develops at the outlet, propagating upwards through a dense shear flow while 
material at the bottom corners of the container remains static. For large enough outlets, the 
discharge flow is continuous; however, by contrast with the clepsydra for which the flow velocity 
depends on the height of fluid left in the container, the discharge rate of granular silos is constant. 
Implementing a plastic rheology in a 2D Navier-Stokes solver (following the /i(I)-rheology or a 
constant friction), we simulate the continuum counterpart of the granular silo. Doing so, we obtain 
a constant flow rate during the discharge and recover the Beverloo scaling independently of the 
initial filling height of the silo. We show that lowering the value of the coefficient of friction leads 
to a transition toward a different behavior, similar to that of a viscous fluid, and where the filling 
height becomes active in the discharge process. The pressure field shows that large enough values 
of the coefficient of friction (~ 0.3) allow for a low-pressure cavity to form above the outlet, and can 
thus explain the Beverloo scaling. In conclusion, the difference between the discharge of a hourglass 
and a clepsydra seems to reside in the existence or not of a plastic yield stress. 

PACS numbers: 45.70.-n, 83.60.Rs, 46.35.+z,83.10.Rs 



I. INTRODUCTION 

Granular matter is well-known for its ability to behave 
like a solid or like a fluid depending on the stress it is sub- 
jected to, transiting from one state to the other over a 
few grain diameters. This phenomenon is best described 
in terms of internal friction: a granular packing can re- 
sist shear stress below the friction threshold and remain 
static, but starts flowing when the friction threshold is 
reached. In this respect, granular matter resembles other 
more classical visco-plastic materials, characterized by a 
yield stress and a viscosity, for instance Bingham plastics. 
The granular silo is one of the many interesting illustra- 
tions of this thixotropic property of granular matter: a 
rapid flow develops at the outlet, propagating upwards 
through a dense shear flow while material at the bottom 
corners of the container remains static [IHS] ■ Beside its 
obvious industrial relevance, the phenomenology of the 
silo discharge is in itself intriguing, and raises an im- 
portant interest from the scientific community. For nar- 
row outlets, arches forms and vanish alternatively, clog- 
ging the flow with a probability depending on the out- 
let dimension, and leading to intermittency [7-10 ; this 
regime, not accessible through continuum modeling, is 
not addressed in this work. For larger outlets, the flow 
is continuous; however, in contrast with the discharge of 
a Newtonian fluid for which the flow velocity depends 
on the height of fluid left in the container, the discharge 
rate of granular silos is constant. The independence of 
the discharge rate on the filling height is best demon- 
strated by the well-known Beverloo scaling, which relates 
the flow rate (Q) to the outlet size (L) [IT]: in a very ro- 
bust manner, experiments and discrete simulations find 
Q = C'y/9(L—kd) N ~ 1 ' 2 , where N is the dimension of the 
problem, d the grains diameter and C and k are constants 



(we consider two-dimensional silos only in the following 
ie N = 2). Because it implies that the velocity of the ma- 
terial flowing from the silo does not scale like the square 
root of the pressure, the Berverloo scaling is commonly 
accepted as the evidence of a screening effect responsi- 
ble for a constant and low value of the pressure (ie lower 
than the hydrostatic prediction) in the area of the outlet 
[12 [T3] . We may suppose that this is why the hourglass 
was used by sailor men for navigation: the oscillations of 
the ship, tilting the device and thereby changing the pres- 
sure inside, probably do not affect the discharge rate as 
much as it would in a clepsydra (or water-clock) , making 
the hourglass more reliable on board. The traditional 
physical explanation for this screening effect resorts to 
Janssen's analysis: the friction forces mobilized at walls 
reduce the apparent weight of the material in the silo 
and prevent the bottom area to sense the pressure, now 
partly sustained by the walls [1 \7\ HHTB]. This effect 
is expected to play a role only if the width of the silo is 
smaller than the filling height. Yet, the Berverloo scaling 
also holds for wide silos |T7j. Moreover, the existence of a 
Janssen screening effect would result in the local pressure 
scaling like the width of the container, yet the Berverloo 
scaling involves only the outlet size. Hence, it seems very 
uncertain that the Janssen effect is responsible for the 
Berverloo scaling and the silo constant discharge rate. 
This was in fact demonstrated in [181 where experiments 
using inclined silos show that the Berverloo scaling holds 
for any degree of tilt up to subvertical values. Experi- 
mental work using horizontal silos also points at a similar 
conclusion [19j . On the other hand, measurements show- 
ing a dip in the value of the pressure in the area of the 
outlet may only reflect the fact that the existence of the 
outlet itself creates a low pressure boundary condition, 
independently of any Janssen screening |13j . This does 
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not mean however that friction is not important during 
the discharge of a granular silo: in this contribution, we 
argue that the role of friction does not involve the walls, 
but the bulk of the flow through the existence of a yield 
stress. 

Implementing a plastic rheology (using either the p(I)- 
rheology [20 or a constant friction) in a 2D Navier-Stokes 
solver [211 122j , we simulate the continuum counterpart of 
the granular silo. Doing so, we observe a constant flow 
rate during the discharge and recover the Beverloo scal- 
ing independently of the initial filling height of the silo. 
However, we show that lowering the value of the fric- 
tion lead to a transition toward a different scaling where 
the filling height becomes active in the discharge pro- 
cess. These results support the idea that the existence of 
a frictional yield stress can alone control the discharge of 
the granular silo without any Janssen effect entering into 
play. 



II. THE CONTINUUM GRANULAR SILO 

The simulations were performed using the Gerris flow 
solver in two dimensions, which solves the Navier-Stokes 
equation for a bi-phasic mixture applying a Volume-Of- 
Fluid approach (5]] [55]. The existence of two fluids trans- 
lates numerically in different properties (viscosity and 
density) on the simulation grid following the advection of 
the volume fraction representing the proportion of each 
fluid. In our case, one fluid stands for granular mat- 
ter (characterized by the coefficient of internal friction 
p) and the other stands for the surrounding air (with a 
lower density and lower viscosity, see [53] for details); the 
position of the interface between the two is solved in the 
course of time based on the spatial distribution of their 
volume fraction. The viscosity rj of the granular matter 
is defined by mean of the friction properties [20] : 

V = min (^J^,Vmax^j , (1) 

where p is the effective coefficient of friction of the gran- 
ular flow, P is the local pressure and D 2 is the second in- 
variant of the strain rate tensor D: D 2 = \fDij Dij ■ For 
large values of D 2 , the viscosity is finite and proportional 
to p and P; when D 2 reaches low values, the viscosity rj 
diverges. Numerically, this divergence is bounded by a 
maximum value ?7 max chosen to be 10 4 times the mini- 
mum value of 77; we have checked that the choice of 77 max 
did not affect the results as long as ^ max is large enough. 
Based on earlier work on continuum modeling of rapid 
non-uniform granular flows showing the better perfor- 
mances of the /x(/)-rheology compared to constant fric- 
tion [53], the effective friction properties p of the gran- 
ular continuum is calculated using the following depen- 
dence: p is a function of the non-dimensional number 
I = dD 2 1 \JPj pi where d is the mean grain diameter and 
p the density (W = 90d in the following, with W the 
silo's width) through the relation p = p 8 + where 
p s , Ap and I are constants J5UM53]. Following [53J, Ap 



FIG. 1: Pressure field during the discharge of a plastic silo of 
width W, normalized outlet size L = 0.125 and normalized 
filling height H = 0.9 (normalized by W) at i = 0, i\ = 0.8, 
I2 = 7.6, £3 = 11.5, and at I4 the final state (normalized by 
\JW I g). The color scale varies from one picture to the other 
to ensure maximum contrast: the highest bound (red color) 
is set to P = 0.6 for to, ti and I2, to P = 0.36 for and to 
P = 0.12 for ii (normalized by pgW). 

and Iq were set respectively to 0.28 and 0.4 and not var- 
ied. The value of the static coefficient of friction p s was 
set to 0.32, 0.2 and 0.1 in order to study its role in the silo 
discharge. Moreover, simulations with a constant friction 
model, ie without dependence on /, were also performed 
(section [v]). 

The silo is flat-bottomed, of width W , filling height H 
and with an outlet of size L (see Figure [TJ. The width 
W is divided in 64 computation cells in the bulk, refined 
to 256 at the bottom, so that the outlet is defined using 
16 to 72 computation cells. In the forthcoming analy- 
sis, all quantities are normalized as follows: silo height 
H = H/W , outlet size L = L/W, volume of material left 
in the s ilo V = V/W 2 , flow rate Q = Q/Wy/gW, time 
t = t/s/W/g and pressure P — P/ pgW . 
A no-slip boundary condition is imposed at the side- walls 
and at the bottom-wall; additional simulations with a 
free-slip boundary condition at the side-walls show that 
the aspects discussed in this paper remain unchanged. A 
zero pressure condition is imposed at the top- wall and at 
the outlet. 



III. A CONSTANT DISCHARGE RATE 

Figure [T] shows the time evolution of a continuum gran- 
ular silo of initial filling height H — 0.9, and outlet size 
L = 0.125; the static friction is set to p s = 0.32 (with 
Ap = 0.28 and Iq — 0.40). The color scale represents the 
pressure field. We observe that the pressure field strongly 
differs from what would be expected in the hydrostatic 
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FIG. 3: Volume V of material remaining in the silo as a 
function of time for L — 0.1875 and filling heights H — 1.5, 
H — 2.5, H = 3.5 and H — 4.5 for a fi(I) plastic flow (fi s = 
0.32, A/x = 0.28 and Io = 0.40). Inset: Berverloo scalings 
corresponding to each case. 
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FIG. 2: a- Volume V of material remaining in the silo as 
a function of time for an outlet size L — 0.125 and filling 
height H = 0.9 in the case of a plastic flow (fj, s — 0.32, A/x = 
0.28 and Io — 0.40) and in the case of a Newtonian flow 
(r) = O.Olpg? W 3 ); b- Berveloo scaling obtained for L varying 
between 0.0625 and 0.28125. 



dh 
h(t) = 



-LV2h, 



r H- —i 
V2 



case, and is non-uniform in the transverse direction. The 
region above the outlet coincides with a low pressure cav- 
ity surrounded by two high-pressure dome-like areas: the 
pressure jump at the outlet overcomes the frictional yield 
stress, creating a large shear and a low viscosity area. 
Meanwhile, pressure gradients decrease in the bulk, and 
a highly viscous mass forms above the outlet. The fact 
that the yield stress is frictional, and depends on the pres- 
sure, implies that it readapts throughout the discharge, 
thereby probably allowing the cavity to survive until the 
end of it. When the material left in the silo stops flowing, 
it remains at equilibrium with a shape depending on the 
yield stress, ie depending on the coefficient of friction p. 
Figure [2]-a shows the volume of material remaining in the 
silo in the course of time for the same system. We observe 
a linear evolution throughout the discharge for the granu- 
lar fluid, revealing a constant flow rate as in real discrete 
granular silos. We measure the value of the viscosity 
in the vicinity of the outlet and find a roughly constant 
value during the discharge equal to rj — 0.01 (normalized 
by pg^W?). Using this value for the viscosity, we sim- 
ulate the discharge of a silo filled with Newtonian fluid; 
the corresponding volume- us-time evolution is shown in 
Figure[2]-a, and is non-linear, suggesting a dependence on 
the height of material left in the silo, as in the case of an 
hydrostatic pressure field. For comparison, we plot the 
solution of the (non-dimensional) Torricelli discharge for 



where h is the instantaneous height of material remain- 
ing in the silo (normalized by W) at time t (normalized 
by yWjg). Torricelli's discharge matches the onset of 
the discharge of the Newtonian fluid provided we chose a 
smaller numerical value for L than that of the simulation 
(L = 0.0714 instead of 0.125). Comparing the discharge 
of the Newtonian fluid and the granular plastic fluid thus 
points at the plastic property of the flow as responsible 
for the constant nature of the discharge rate in the latter 
case. 

We observe that the flow rate remains constant through- 
out the discharge of the plastic silo over a large range of 
outlet size L. Varying L between 0.0625 (ie 16 compu- 
tation cells) and 0.28125 (ie 72 computation cells), we 
measure the flow rate Q, and search for a relation satis- 
fying the shape of the Beverloo scaling: 

Q = C(L-ky (2) 

where C and k are constants. For a continuum silo, the 
numerical value of k is expected to be zero, in contrast 
to granular silos where the grain diameter imposes a vol- 
ume of exclusion reducing the effective size of the outlet. 
Imposing k = 0, we recover the Berverloo scaling with 
a good accuracy, giving C = 1.4 (see Figure [2j-b). Note 
however that making no assumption on the fitting pa- 
rameters, the best fit gives k = 0.00938 = 0.85J, where 
d = 1/90 is the grain diameter used in the /i(/)-rheology. 
Although the value of k is not completely negligible, we 
consider nevertheless that it is not physically significant. 



4 




Ms = 0.32 , A|J = 0.28 Ms = 0.32 , Am = 0.28 



FIG. 4: Pressure field in the early stage of the discharge for 
two granular plastic silos of initial filling heights H = 1.5 
(right) and H = 4.5 (left), outlet L = 0.1875 and /-dependent 
friction coefficient fx s = 0.32 (A/i = 0.28 and I = 0.40). The 
color scale is identical on both images: the highest bound 
(red color) is set to P = 1.3; the pressure jump between two 
isolines is 0.15. 



IV. INCREASING THE FILLING HEIGHT 

To check whether the discharge rate remains constant 
irrespective of the initial filling height, we perform 
additional simulations with H = 1.5, 2.5, 3.5 and 4.5, 
with the same rheological parameters as previously 
(jj,„ = 0.32, A/i = 0.28 and J = 0.40). Figure § shows 
the volume V left in the silo in the course of time for an 
outlet size L = 0.1875 for all four cases. We observe a 
linear evolution of very similar slope, suggesting that the 
discharge rate is essentially constant and independent 
of the filling height. However, closer inspection shows 
that the slope varies slightly during the discharge: for 
H = 4.5 for instance, the initial flow rate has decreased 
of 1.4% halfway through the discharge. Moreover, a 
slight increase of the flow rate is observed for larger 
filling height: approximating the discharge by an affine 
function over its full duration, we find a flow rate 
increase of 1.8% for H = 2.5 and an increase of 3.9% for 
H = 4.5 compared to the case of H = 1.5. 
Measuring the flow rate Q in the early stage of the 
discharge for the different values of H and different 
outlet dimension Z, we recover the Berveloo scaling 
but with coefficients varying slightly with the value of H 
(Figure [3j inset). This weak influence of the initial filling 
height is maximum in the early stage of the discharge, 
but vanishes at the end: the curves shown in Figure 
[3] can be superimposed when shifted towards the final 
stage of the discharge. 

Altogether, the flow rate during the discharge of con- 
tinuum granular material is very weakly affected by the 
initial filling height H . By contrast, the pressure field 
strongly changes with the value of H, as is visible on 
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FIG. 5: Volume of material V remaining in the silo as a 
function of time i for L — 0.1875 and filling height H = 4.5, 
for a plastic flow with fi a = 0.10, fi a = 0.20, and fi s = 
0.32 (A/i = 0.28 and I = 0.40). Inset: Flow rate Q as a 
function of outlet size L for /t 3 — 0.10 and fi s — 0.32. 
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FIG. 6: Pressure field in the early stage of the discharge 
for a granular plastic flow with a constant friction coefficient 
(is = 0.32 (A/i = 0., left), a /-dependent friction /i s = 0.10 
(A/i = 0.28 and /o = 0.40, center) and for a Newtonian flow 
of viscosity fj — 0.01 (normalized by pg^Ws) (left). Outlet 
size L — 0.1875, filling height H — 1.5. The color scale is 
identical on all images and identical to Figure [4] the high- 
est bound (red color) is set to P = 1.3; the pressure jump 
between two isolines is 0.15. 

the two snapshots shown in Figure [4] for H = 1.5 and 
H = 4.5 in the early stage of the discharge. In both 
cases however, we observe a low pressure cavity in the 
vicinity of the outlet, suggesting that the discharge is 
affected by this local pressure condition and insensitive 
to the mean pressure in the silo. 



V. INFLUENCE OF THE INTERNAL 
FRICTION 

We suspect the deviation of the discharge rate from the 
hydrostatic case to be the result of the non-newtonian na- 
ture of the ^(I)-rheology, and more specifically of the ex- 
istence of a yield stress; accordingly, we expect the value 
of the coefficient of friction to have important repercus- 
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sions on the discharge flow. Figure [5] shows the discharge 
of a silo with H = 4.5 and outlet L = 0.1875 for differ- 
ent values of the coefficient of static friction fi a — 0.10, 
(a s = 0.20 and /i s — 0.32: we observe that a smaller 
friction coefficient induces an earlier departure from the 
linear evolution observed for larger friction. Since fric- 
tion does not only act upon the yield stress, but also 
upon the value of the viscosity, it affects the discharge 
duration too. 

For the three values of the static coefficient of friction, 
and for a silo with H = f.5, we vary the outlet size L 
between 0.0625 and 0.28f25 and measure the flow rate 
Q. For fi s — 0.10 and fi s — 0.20, the flow rate is not 
constant all through the discharge. Hence, we consider 
the early stage of the discharge only, where an affine re- 
gression seems reasonable. As already seen in previous 
sections, the flow rate in the case (i s — 0.32 obeys the 
Beverloo scaling (see Figure [2]) ; in the case of lower fric- 
tion however (fi s = O.fO), the Beverloo scaling is no more 
relevant, and the flow rate increases linearly with the out- 
let size L (Figure [5J inset), thus making a dependence on 
H 1 / 2 dimensionally possible. 

An other aspect is the importance of the / dependence on 
the friction model. Although fully addressing this ques- 
tion implies detailed comparison with discrete systems 
(either experimental or numerical), we can nevertheless 
check how the height independence of the silo discharge is 
affected by suppressing this dependence and set Afi = 0. 
Figure [6] shows the pressure field in a silo with fi s = 0.32, 
A/i = and H — 1.5 (to compare with Figure we 
observe no significant difference with the case Afi ^ 0. 
Plotting the discharge rate against outlet size allows for 
the recovery of the Beverloo scaling as well (not shown). 
It is thus clear that the existence of a frictional threshold 
alone is enough for reproducing the height independence, 
without the need of a / dependence. The effect of the 
latter being to increase the friction close the the outlet, 
it may become important when comparing quantitatively 
discrete and continuum granular silo discharge rate. This 
discussion is however beyond the scope of this paper. 
Figure [6] also presents the pressure field in a plastic silo 
with weak friction fj, s =0.1 and in a Newtonian silo of 
viscosity fj — 0.0 1 (normalized by pg^W~s), and shows 
that the existence of the low pressure cavity is contingent 
on the existence of a yield stress induced by sufficiently 
large friction (Figure [6]). 

VI. DISCUSSION 

Granular matter forms a specific class of plastic fluids 
for which yield stress and viscosity are not independent, 



but related through friction properties. Hence discrim- 
inating between the respective roles of yield stress and 
viscosity is difficult. However, the continuum simulations 
presented in this paper show that i) a friction-dependent 
viscosity (either fi(I) or constant friction) leads to a con- 
stant flow rate as is the case for real granular silo and 
ii) decreasing the coefficient of friction leads to a depen- 
dence on the height of material remaining in the silo, as 
is the case for Newtonian fluids. Decreasing the coeffi- 
cient of friction to a value as small as 0.1 may not be 
physical for granular matter; however, this extrapolation 
tends to show that the constant shear rate results from 
the existence of a yield stress. For sufficiently large val- 
ues of the coefficient of friction (ie > 0.3), the pressure 
field shows the existence of a low pressure cavity in the 
vicinity of the outlet, despite high values of the pressure 
elsewhere in the silo (Figure [4]) . The fact that the yield 
stress (^t-P) depends on the mean pressure in the system 
is certainly crucial in the ability of the cavity to remain 
identical all though the discharge. This cavity does not 
exist during the discharge of a Newtonian flow, and is 
very reduced in the case of a plastic flow with a small 
friction coefficient fi s = O.fO (Figure [6]). It thus seems 
reasonable to conclude that the frictional properties of 
the continuum granular material are at the origin of the 
low-pressure cavity forming above the outlet; this coin- 
cides with a high-shear-low-viscosity region resembling 
the free fall arch described in granular systems [T^] . The 
fact that the flow is controlled by the local conditions 
close to the outlet rules out the Janssen effect as expla- 
nation for the Beverloo scaling (as did before [TBI US] 
experimentally). In conclusion, the difference between 
the discharge of an hourglass and a clepsydra seems to 
reside in the existence or not of a frictional yield stress. 
The ability of the continuum /i(/)-rheology to reproduce 
the behavior of granular systems when implemented in a 
Navier-Stokes solver was discussed in [23]. In this con- 
tribution, we show the ability of the same model to re- 
produce the silo phenomenology, without discussing in 
depth the relevance of the /-dependence. It should be 
noted that a constant friction model leads to the silo 
phenomenology too, as far as height independence and 
Beverloo scaling general shape are concerned. More de- 
tailed aspects like comparison of velocity and pressure 
profiles, inner deformations and surface shape are beyond 
the scope of the present work. They are the subject of a 
forthcoming dedicated communication. 
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